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Abstract. We report the results of self-gravitating simu- 
lations of spiral galaxies, modeled by stellar and gaseous 
components, developed to investigate in particular the role 
of dissipation in the evolution of galaxy disks. The gas 
disk is simulated by the Beam-Scheme method, where it 
is considered as a self-gravitating fluid. The results suggest 
that the gravitational coupling between the stars and gas 
plays a fundamental role in the formation and dissolution 
of stellar bars, depending on the gaseous mass concentra- 
tion and on the degree of dissipation. In addition we re- 
mark that initially concentrated gas disks can be unstable 
to the one-armed (m=l) spiral perturbations, which may 
explain the lopsided features observed in the gas distribu- 
tion of the late-type isolated galaxies. The development of 
the m=l feature slows down the radial gas flows towards 
the center, since the large-scale gravity torques are then 
much weaker. 

Key words: Galaxies: evolution; Galaxies: general; 
Galaxies: ISM; Galaxies: kinematics & dynamics; Galax- 
ies: spiral; Galaxies: structure 



1. Introduction 

While the gravitational instabilities of a purely stellar disk 
have been widely studied, mainly through simulations, 
and begin to be well understood, the situation is much 
more open for a self-gravitating disk composed of stars and 
gas. Even when a stellar disk is stabilised against axisym- 
metric instabilities (Toomre criterion, 1964), it can be the 
site of spiral and bar instability (e.g. Sellwood & Wilkin- 
son 1993), together with z-instabilities (e.g. Combes et al 
1990). These heat considerably the stellar disk, that can 
then no longer sustain spiral structure. Instabilities can 
be suppressed by reducing the effective self-gravity of the 
stellar disk, either through the addition of a hot bulge 



or halo (Ostriker & Peebles 1973), or through disk heat- 
ing, i.e. increasing the velocity dispersion (Athanassoula 
& Sellwood 1986). 

The behaviour of a self-gravitating disk of both stars 
and gas, even if the latter represents only a few percents 
of the total mass, presents much more variety and com- 
plexity. Due to the dissipation, the gaseous component 
remains cool, and evacuates the heating due to the gravi- 
tational instabilities: spiral structure can then be contin- 
uously renewed (e.g. Sellwood & Carlberg 1984), even in 
the stellar disk. Criteria for stability are more complex, 
the gravitational coupling between gas and stars making 
the ensemble unstable, even when each component would 
have been separately stable (e.g. Jog & Solomon 1984, 
Bcrtin & Romeo 1988, Elmegreen 1994). Gravity torques 
exerted on the gas produce strong radial flows, and gas can 
accumulate at Lindblad resonances and form several rings 
(e.g. Schwarz 1981, Combes 1988). The strong gas con- 
centration in the nucleus can destroy stellar bars (Friedli 
& Benz 1993), after having sometimes triggered a nuclear 
bar within the main one (Friedli & Martinet 1993). This 
bar destruction can be explained by an increase of chaotic 
orbits in the stellar bar (Hasan & Norman 1990). The co- 
pious inflow of gas into the nuclear region may be one of 
the main source of nuclear activity in galaxies (Phinney, 
1994). 

When the gas possess too much self-gravity, however, 
it forms lumps through Jeans instability, and the lumpi- 
ness of the gas can scatter the stars, randomize their 
motions, and prevent any bar formation (Shlosman & 
Noguchi 1993). While the gas is triggering bar instability 
when it represents only a few percents of the total mass, 
it can play the inverse role when its mass is above « 10%. 

All the complexity of the evolution of a self-gravitating 
galaxy disk composed of stars and gas has not yet been 
explored; moreover, the detailed physics of the interstellar 
medium has not yet be quantified: we still do not know 
the exact amount of dissipation, or the energy balance 
between dissipation and star formation, etc.. We want to 
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study to what extent the galaxy evolution depends on 
these parameters, through the gas behaviour. We want 
to compare the behaviour of the two main gas phases, 
condensed and diffuse. 

Observations show that the gas component is a multi- 
phase medium, spanning a wide range of densities and 
temperatures. To drastically simplify, we can consider two 
distinct forms - dense, cold and clumpy molecular clouds 
and warm diffuse gas - with different spatial distribution 
but having comparable total masses. The two aspects of 
the gas component have been studied in the literature: en- 
sembles of interstellar clouds that undergo collisions, but 
do not behave as a fluid (Schwarz 1981, 1984; Roberts & 
Hausman 1984; Combes & Gcrin 1985). The collisions take 
place mainly along spiral arms. This method describes 
essentially the molecular gas component. The response 
of the diffuse gas in barred potentials was also explored 
(Huntley et al 1978; Roberts et al. 1979; Sanders & Tubbs 
1980; van Albada & Sanders 1982; Prendergast 1983; van 
Albada 1983; Contopoulos et al. 1989). In these calcula- 
tions the gas was considered as a fluid, submitted to pres- 
sure forces, and undergoing shocks. In most cases, the gas 
self-gravity was not considered. 

In the present paper we concentrate on the self- 
gravitating diffuse phase, and present numerical simula- 
tions of spiral galaxies modeled by two components, stars 
and gas-fluid. The primary motivation is to quantify the 
amount of dissipation (cither artificial, turbulent or "grav- 
itational" ) and radial flow of the gas in the formation and 
evolution of the gravitational instabilities observed in spi- 
ral disks. In particular we discuss phenomena that can 
be explained by the star-gas coupling. Future studies will 
compare the condensed with the diffuse phase behaviour, 
and will study a combination of the two. 

This paper is organized as follows. In § 2 we discuss 
the computational methods, the galaxy model used and 
the initial conditions of our simulations. In § 3 we present 
the results concerning the m=2 instabilities, and in § 4 the 
m=l and one-armed spiral instabilities. Section 5 quanti- 
fies the effects of viscosity and section 6 summarizes our 
conclusions. 

2. Numerical Simulations 

2.1. Computational Methods 

We simulate the self-gravity of both the stellar and gaseous 
components, and solve the N-body problem by the 2D- 
Fast Fourier Transforms method. At each integration step, 
density (of stars and gas) is computed at each point of a 
grid of 256 x 256 square cells. Fourier transform of density 
and inverse transform of potential are then carried out. 
The potential is softened within 250 pc, which is also the 
cell size. Each particle contributes to the density of the 
four nearest grid points, according to cloud-in-cell (CIC) 
interpolation procedure (Birdsall & Fuss 1969). There are 



36000 particles representing the stellar component. The 
gaseous component is represented by a fluid. 

The fluid is simulated according to the Beam-Scheme 
method developed by Sanders & Prendergast (1974). This 
technique for solving the hydrodynamical equations, is 
an explicit Eulerian scheme. The gas is evolved on a 
256 x 256 2D-Cartesian grid, with mass, momentum, and 
energy densities being specified in each cell. The presumed 
Maxwellian velocity distribution in each gas cell is approx- 
imated by five delta functions (or beams), which reproduce 
the appropriate moments of the velocity distribution func- 
tion. These beams are allowed to move over a time step 
transporting mass, momentum and energy into adjacent 
gas cells. Transport is taken into account to determine the 
new mass, momentum, and energy in each cell, and these 
moments are used to describe the new Gaussian velocity 
distributions for each cell of the gas. 

2.2. Galaxy Model 

As for the initial conditions, the model galaxy is com- 
posed of a stellar disk, an analytical spherical bulge and a 
gaseous disk. The stellar particles are initially distributed 
according to a Toomre disk (Toomre, 1962) of mass M<j 
and radial scale length a^. The gaseous component of mass 
M g is modeled by another Toomre disk, of radial scale 
length a g , with velocities distributed in beams according 
to the Beam-Scheme method. A Plummer potential (Bin- 
ney & Tremaine 1987) of radial scale length at, represents 
the bulge. In all simulations the bulge and the disk have 
the same mass (Md = Mb) and the gas mass corresponds 
to 10% of the total mass of the galaxy (M g = 0.1M t ; 
M t = 10 11 M ). The initial cutoff radius of the stellar 
and gaseous disks are 12 kpc and 16 kpc, respectively. 
The characteristic scale of the stellar disk, of the bulge 
and of the gaseous disk are free parameters. Their values 
are shown in Table 1. Initially, the stellar particles have a 
velocity dispersion as a function of radius corresponding 
to the Toomre parameter Q = 1 all over the disk. The 
rotational velocity is corrected for asymmetric drift. The 
units of time, distance and mass for all simulations are 
10 7 years, 1 kpc and 4 x 10 9 Mq respectively. All the ex- 
periments have been evolved through a total time equal 
to 2.2 x 10 9 years. The time-step is 0.1 unit= 1 Myrs. 

In the experiments A, B, C and D the gaseous compo- 
nent is not included. These are control simulations in order 
to probe the gas influence. In the E, F, G and H exper- 
iments the stellar components arc described by the same 
initial conditions as in the A, B, C and D experiments 
but a gaseous disk, of scale length a g — 6 kpc, is included. 
A comparative analysis of these experiments allow us to 
draw important conclusions about star-gas coupling. Ex- 
periments I and J test the role of the initial concentration 
of the gaseous disk in the evolution of the system. 
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Table 1. Parameters of the Experiments 



Exp. 


a-d 


a b 


a g 


A 


3.0 


1.0 




B 


4.5 


1.0 




C 


3.0 


2.0 




D 


4.5 


2.0 




E 


3.0 


1.0 


6.0 


F 


4.5 


1.0 


6.0 


G 


3.0 


2.0 


6.0 


H 


4.5 


2.0 


6.0 


I 


3.0 


1.0 


4.0 


J 


3.0 


1.0 


8.0 



3. M=2 instabilities 

In this section we discuss the role of gas dissipation in the 
formation and maintenance of gravitational instabilities 
(m = 2 spirals, bars). 

3.1. Influence of the gas 

In figure [j] we present the linear isodensity contours of 
the stars for the experiments A, B, C and D at the time 
t=140 (10 7 yrs). In these four experiences, as described 
in the preceding section, the gaseous component is not in- 
cluded. We can note that in spite of the large bulge mass 
(Md = Mf,) a central bar develops. The spiral arms formed 
initially are only transient. They help to transfer angular 
momentum outwards to make the bar grow. The charac- 
teristic time-scale for bar formation in these experiments 
varies according to the degree of concentration of the ini- 
tial bulge and stellar disk. Models with a smaller bulge 
scale length at,, developped their bar in a shorter time- 
scale than those with a larger initial value of ab- This can 
be checked by comparing experiments A and C or the ex- 
periments B and D, and can be interpreted by a shorter 
dynamical time-scale in the center of A and C. The bar is 
also shorter in length in Exp. A and C, and the corotation 
radius is smaller. Considering now the same initial mass 
concentration for the bulge (experiments A and B) we ver- 
ify that disks initially more concentrated (experiment A) 
allow a more rapid bar formation. This could be due to the 
larger self-gravity of the stellar disk in the center, together 
with the shorter td there. 

The main conclusion from these four experiments is the 
high longevity of the bar: the bar kept its amplitude up to 
the end of simulations (2.2 Gyrs). In the absence of gas, 
bars are long-lived, even in the presence of a concentrated 
bulge. 

On the contrary, a great variety of scenarios can occur 
when the gaseous component is included. In figures |^, [|, 
H and |^ we present the linear isodensity contours of the 
stars and gas for the experiments E, F, G and H, at four 
different epochs (t=20, 80, 140 and 200 x 10 7 yrs). In 



3 

these experiments, the stellar and gaseous components can 
develop a bar (experiment E), in the same way as pure 
stellar runs, or evolve to various degrees of central mass 
concentration (Exp. F, G and H). In that case either the 
bar forms, and is then destroyed rapidly (experiments G 
and H) or even the bar does not form at all (experiment 
F). 



Fig. 1. Exp. A, B, C and D: linear isodensity contours of the 
stars at epoch t=140 (in 10 7 yrs). X and Y are in kpc. 



Fig. 2. Exp. E: linear isodensity contours of the stars (left 
side) and of the gas (right side) at epochs t=20, 80, 140 and 
200 (in 10 7 yrs). X and Y are in kpc. 
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Fig. 3. Exp. F: linear isodensity contours of the stars (left 
side) and of the gas (right side) at epochs t—20, 80, 140 and 
200 (in 10 7 yrs). X and Y are in kpc. 



Fig. 5. Exp. H: linear isodensity contours of the stars (left 
side) and of the gas (right side) at epochs t=20, 80, 140 and 
200 (in 10 7 yrs). X and Y are in kpc. 



Fig. 4. Exp. G: linear isodensity contours of the stars (left 
side) and of the gas (right side) at epochs t=20, 80, 140 and 
200 (in 10 7 yrs). X and Y are in kpc. 



3. 2. Bar amplitude and pattern speed 

As described above, the bar strength and evolution de- 
pends on the central mass concentration and the gas frac- 
tion. We focus in this section on the bar evolution, the 
results of the experiment F where a central mass concen- 



tration is obtained without the formation, even transient, 
of a bar will be discussed in the next section. 

To better identify the m = 2 perturbation, we have 
computed, for all simulations, the angular velocity of the 
pattern Q p as a function of radius. A mode would be 
characterized by a constant £l p over a range of radii. A 
Fourier analysis of the density distribution or of the po- 
tential allows to decompose the perturbation present in 
the galaxy in several components (m). The index m in- 
dicates the symmetry of each component or, in the case 
of spirals, the number of arms. The analysis provides, at 
each epoch, both the strength or amplitude of the m com- 
ponent, which growth rate can then be estimated and its 
phase, from which Q p (r) can be derived. We have com- 
puted the growth rate for the spiral-bars in the exper- 
iments A, E, G and H assuming exponentially growing 
modes. The behavior of the bar amplitude, the maximum 
ratio of the m = 2 tangential force to the radial force (P2) 
as a function of time is shown in figure |. While for the 
experiments A and E p2 remains approximately constant 
as in a steady-state (although the bar weakens somewhat 
in experiment E), in experiments G and H it decreases 
rapidly. The growth time-scales obtained for the bars are 
approximately 7 t d (Exp. A), 6 t d (Exp. E), 5 td (Exp. G) 
and 4 td (Exp. H) where t d is the dynamical time scale 
(« 4x 10 7 yrs). We remark that in Exp. A and E the 
behavior of pi in the first instants of the simulation (be- 
fore the exponential bar growth) indicates the presence 
of transient spiral structure preceding bar formation. We 
have calculated for the experiments A and E the angular 
velocity (or pattern speed) for the bar perturba- 

tion. In figure [7| are shown the f2 (solid curve), £1 + k/2 
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and f2 — k/2 (dashed curves) frequency curves, where n 
represents the epicyclic frequency and fl is the angular 
velocity of the galaxy. The expected position of the Lind- 
blad resonances are indicated at the top. It can then be 
noted that the bars are confined within their corotation 
radius in both experiments. 



Fig. 6. The strength of the bar perturbation p2 as a function 
of time (in 10 7 yrs units) for the experiments A, E, G and H. 
The m=2 amplitude has been averaged over radius to 5kpc. 



Fig. 7. The frequency curves (in units of 100 km/s/kpc) as 
a function of radius (in kpc) and the position of the pattern 
speed Q p (r) for experiments A (top) and E (bottom). 



Fig. 8. Exp. E, F, G and H : surface density (in code units) 
of the stars and g function of radius (in kpc) at four 

different epochs t=20, 80, 140 and 200 (in 10 7 yrs). The letters 
at the bottom left corner of each figure indicate (s) stars and 
(g) gas. 



Figure || displays the evolution of the surface density of 
stars (left side) and gas (right side), as a function of radius, 
at four different epochs. These quantify the consequences 
of radial gas flows and mass concentration of stars and gas 
in the central region. Radial flows correspond to angular 
momentum transfer driven by non-axisymmetric gravita- 
tional instabilities, and the subsequent central mass con- 
centrations react back on gravitational instabilities. 

As is apparent in figure |§|, in all simulations the sur- 
face density of stars still follows, at the first epoch (t=20), 
the imposed initial conditions. On the contrary, the gas 
surface density is determined at the first epoch, by the 
bulge mass concentration. Such behavior determines the 
future evolution of the stellar density distribution. In other 
words, at t—2Q the gas, in the experiments E and F 
(at, = 1), is more concentrated than in the experiments 
G and H. This strong initial concentration in experiment 
E is the cause of a high central accumulation of stars, and 
after some time favors the formation and maintenance of 
the bar. More precisely, it avoids a violently unstable so- 
lution, that will heat the disk and destroy the bar. In 
experiments G and H the bar forms more rapidly but it is 
not maintained. At t=80, 140 and 200 we observe that the 
gas in experiments E, G and H can reach the same max- 
imum of density. While in experiment E the region out 
of the bar is depopulated from its gas, the contrary oc- 
curs in experiment F. In this experiment a strong central 
concentration of gas at the initial epochs together with a 
considerably lower concentration of stars induce the for- 
mation of an m=l spiral mode in the central region of the 
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galaxy. The concentration of gas and stars is then almost 
completely stopped. 



agreement between the spiral arms developed in the gas 
and in the stars. At t=80, 140 and 200 there is in the gas 
a phase-shift between the surface density and potential, 
implying a large torque felt by the gas for a large range of 
the radii; after corotation, the stellar surface density leads 
the potential, while just before CR the potential leads the 
two components. Corotation is where the potential phase 
shows the maximum gradient. At t=200, a gas spiral fea- 
ture just outside corotation introduces important phase 
shifts, but they are not related to the bar, as shown by di 
the m=2 amplitude in the surface density. 



Fig. 9. Exp. E : the phase (in radian) of m=2 Fourier com- 
ponent for the stellar (left side) and the gaseous (right side) 
surface density and the total potential as a function of radius 
(in kpc) at four different epochs. Also shown is d,2, the ampli- 
tude of the m=2 component of the surface density of stars and 
gas. The scale of d,2 is in the right side. 

We can observe that although the bar obtained in Exp. 
A shows approximately the same growth rate as that ob- 
tained in Exp. E, it is relatively longer. This is certainly 
due to the larger mass concentration of Exp. E, due to the 
gas dissipation. Also spiral structure is present in Exp. E 
in contrast to A, in the gas as well as in the stars. We re- 
mark also that there is a phase shift between the gas bar 
(leading) and that of the stars of approximately 30°in the 
first epochs. This explain the strong gravity torques ex- 
erted by the potential on the gas, which is driven inwards 
(cf Combes 1988). In her study about the secular evolution 
of spiral galaxies, Zhang (1995), suggested the existence of 
a collective dissipation mechanism responsible for the sec- 
ular evolution of the stellar disks of spiral galaxies. Since 
the disk is not in complete equilibrium, there exists also a 
phase shift between the self-consistent spiral potential and 
the stellar density. This phase shift indicates that there 
is a torque applied by the spiral potential on the stars, 
and a secular transfer of energy and angular momentum 
between the stellar disk and the spiral density wave. In 
figure |p we show, for experiment E, the phase of the m=2 
Fourier component for the gas and stellar surface densities 
and for the total potential as a function of radius at four 
different epochs, t=20, 80, 140 and 200 (in 10 7 yrs). At 
the beginning (t=20) gaseous and stellar surface densities 
are in phase with the potential. See also in hgure | the 



Fig. 10. Exp. I, E, and J: linear isodensity contours of stars 
and gas at epoch t=140 (in 10 7 yrs). X and Y are in kpc. 

In order to probe the conditions required for the forma- 
tion of a bar (as in experiment E) or an m=l spiral mode 
(Exp. F), two other experiments (I and J) were done. In 
these experiments we used the same initial conditions for 
the stellar disk as in experiment E but the gaseous com- 
ponent has different degrees of concentration determined 
by the characteristic radius a g . In figure |3~0| we compare 
the linear isodensity contours of the stars and gas at the 
time t=140 for Exp. I, E and J: when the gaseous disk 
is initially less concentrated a bar can form (Exp. J and 
E) while when it is initially more concentrated (Exp. I, 
a g = 4) an m=l spiral mode develops. In this latter case 
the stars and gas concentrations show a final behavior very 
different from that observed in the bar experiments. The 
results can be interpreted comparing the surface density 
of stars and gas at different epochs as shown in figure pA\ 
Although the stellar density of Exp. I does not stand out 
from other experiments at any time, the gas density is 
different. While in Exp. E and J the bar forms, the con- 
centration of gas increases and the region out of the bar is 
depopulated, in Exp. I, the development of an m=l spiral 
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In this section we concentrate on perturbations with 
azimuthal wavenumbcr m=l. The results of one of these 
simulations (Table 1, experiment F), which shows gravita- 
tional instabilities of this nature, is presented in figure [| 
The linear isodensity contours of the gas show features in 
the central region that can be identified as an m=l spi- 
ral mode. The Fourier analysis allows us to reconstruct 
the density distribution of the gas from the combination 
of the various components and single out the contribu- 
tion of each one in the final configuration. In figure [jj] we 
present the linear isodensity contours of the gas distribu- 
tion reconstructed from the two main components (m=l 
and m=2). Comparing figures || and |l2| we can check the 
relative contribution of these components and how they 
add to each other. 



Fig. 11. Exp. E, I, and J : surface density (in code units) of 
stars and gas as a function of radius (kpc) at four different 
epochs, t=20, 80, 140 and 200 (in 10 7 yrs). The letters at the 
bottom left corner of each figure indicate (s) stars and (g) gas. 



mode slows down the gas concentration. We discuss now 
in more details the formation of these m=l spiral modes. 



4. M=l instabilities 

Eccentric asymmetries in the distribution of light in spiral 
galaxies have been known for a long time (Baldwin et al. 
1980; Wilson & Baldwin 1985). Lopsided features are pref- 
erentially observed in the distribution of gas in late-type 
spiral galaxies. In several cases these features can be iden- 
tified as one-armed spirals (m=l mode). More frequently, 
nuclei of galaxies are observed displaced with respect to 
the gravity center, as in M33 and M101 (de Vaucouleurs & 
Freeman 1970). Barred spiral galaxies can have their kine- 
matical center displaced from the bar center (Christiansen 
& Jefferys 1976, Marcellin & Athanassoula 1982, Duval & 
Monnet 1985). The nucleus of M31 reveals such an off- 
centring (Bacon et al 1994) which has been interpreted 
in terms of an m=l perturbation (Tremaine 1995). Miller 
& Smith (1992) have studied through N-body simulations 
of disk galaxies, a peculiar oscillatory motion of the nu- 
cleus with respect to the rest of the axisymmetric galaxy. 
They interpret the phenomenon as an m=l instability, a 
density wave in orbital motion around the center of mass 
of the galaxy. Weinberg (1994) shows that a stellar sys- 
tem can sustain weakly damped m=l mode for hundreds 
of crossing times. A fly-by encounter could excite such a 
mode, and explain off-centring in most spiral galaxies, or 
the forcing of a halo could maintain a time-dependent po- 
tential in the galaxy (Louis & Gerhard 1988). 



Fig. 12. Exp. F: linear isodensity contours do gas recon- 
structed from the m=l and m=2 Fourier components epochs 
t=20, 80, 140 and 200 (in 10 7 yrs). X and Y are in kpc. 

For this experiment, we have calculated the angular ve- 
locity f2 p (r) of the m=l wave. The m=l wave is confined 
in the central region (r < 3 kpc), with a pattern speed of 
the order of 400km/s/kpc. Figure |l3| displays the f2 (solid 
curve), fl + k and fl — k (dashed curves) frequency curves, 
together with the expected position of the resonances (CR 
and OLR). The m=l perturbation develops between the 
center of the galaxy and its OLR radius, slightly above 
3 kpc. Outside this radius, independent perturbations, of 
symmetries m=l and 2, develop (see figure |l4]. The pat- 
tern speed of these is of the order of 10-20 km/s/kpe. In 
figure [l4| we show the linear isodensity contours of the gas 
in Exp. F at twelve consecutive epochs (from t=140 to 
142.75, in 10 7 yrs). We can note that the central region 
(about 3 kpc around the densest cell) of the m=l spi- 
ral wave rotates rapidly, completing almost two rotations 
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Fig. 13. Exp. F: The angular velocity £l p (r) of the m=l wave 
and the frequency curves (in units of 100 km/s/kpc) as a func- 
tion of radius (in kpc). 



Fig. 14. Exp. F: linear isodensity contours of the gas at 12 
close- by epochs (in 10 7 yrs). The circles show the region of 3 
kpc around densest cell. The axis X and Y are in kpc. 



during this time interval while the outer spiral structure 
evolves more slowly. The value of the central fi p (r), consid- 
erably greater than in the outer regions of the galaxy, pre- 
cludes any non-linear coupling of the patterns in the two 
regions. In figure [l^ we show the same linear isodensity 
contours displayed in figure |lj rotated with a constant an- 
gular velocity (S!p(r)=400km/s/kpc). The position of the 
m=l wave at each time in the central region confirms that 
the pattern speed is in fact considerably greater there. 

It is important to remark that, at the initial instants 
of the simulation (t=20) in the central region (r=3 kpc), 
the mass of the gas increases to approximately 6 times its 
initial value, just before the appearance of the m=l mode. 
After that, the growth rate of gas mass in this region slows 
down steadily with time. In the other experiments where 
the m=l does not occur, the gas mass continues to grow 
rapidly. The slowing down of radial gas flow coincides also 
with a stabilisation of the gas angular momentum in the 
center, as soon as the m=l arm is identified. 

Adams et al (1989) have discovered numerically sim- 
ilar gravitational instabilities in gaseous disks associated 
with young stellar objects: eccentric m=l modes, when 
the star did not lie at the center of mass of the system. 
Shu et al. (1990) presented an analytical description of a 
modal mechanism, the SLING amplification, or "Stimula- 
tion by the Long-range Interaction of Newtonian Gravity" . 
This mechanism uses the corotation amplifier, where the 
birth of positive energy waves outside strengthens neg- 
ative energy waves inside. In the SLING mechanism, a 
feedback cycle is provided by four waves outside corota- 
tion; long-trailing waves propagate from the OLR inward 
to the Q-barrier at CR where they refract in short trail- 



Fig. 15. The same linear isodensity contours of the gas at 
12 close-by epochs (in 10 7 yrs) shown in figure 14 rotated of 
D. p (r) = 400km/s/kpc. 

ing waves. These propagate outwards, cross the OLR, and 
reflect back at the outer edge of the disk. This is a criti- 
cal point, the whole amplification mechanism depends on 
the reflecting character of this outer edge. Then a short 
leading wave propagates inwards from the edge, through 
OLR, towards CR, where it refracts again into a long lead- 
ing wave, that is then reflected at OLR. An essential point 
here is also the ability of waves not to be absorbed at OLR, 
that is why this mechanism apply to gaseous disks only. 
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Fig. 16. Exp. F: position of the center of mass of the gas 
(circles) and the center of mass of the stars (*) relative to 
the center of mass of the system (crosses) calculated within 
radii from r=l kpc to 10 kpc at different epochs indicated at 
the top left of each figure (in 10 7 yrs). X and Y are in kpc. 
The position of the center of mass of the gas and of the stars 
calculated within the same radius are connected by a straight 
line. 



Using a WKBJ analysis, they derived from the dispersion 
relation, and the condition of a constructive reflection, the 
required pattern speed for these modes. The instability 
arises uniquely from the displacement of the central star 
from the center of mass of the system, which creates an ef- 
fective forcing potential. In our simulation the same effect 
is occurring, i.e. the center of mass of the gas and of the 
stars is displaced from the center of mass of the system 
and they are displace in opposite positions, inducing the 
m=l wave formation. In other words, the modal mecha- 
nism proposed by Shu et al. (1990) may be in action in 
our experiment. In figures [l6] and 17 we show the position 
of the center of mass of stars and gas calculated within 
radii from r=l kpc to r=10 kpc, at 12 different epochs 
(indicated at the top left of each frame) for experiments F 
and E respectively. In this figures the circles indicate the 
position of the center mass of the gas, the symbol (*) the 
position of the center mass of the stars and the crosses the 
center mass of the system. From the comparison of figures 
|l6| and 17, we note that when there is a bar (Exp. E) the 
position of the center of mass of the gas and of the stars 
agree with the position of the center mass of the system, 
while when an m=l spiral wave is identified, the center of 
mass of the gas and of the stars are in opposite locations 
relative to the center of mass of the system. 



Fig. 17. Exp. E: position of the center of mass of the gas 
(circles) and the center mass of the stars (*) relative to the 
center of mass of the system (crosses) calculated inside radii 
from r=l kpc to 10 kpc at different epochs indicated at the 
top left of each frame (in 10 7 yrs). X and Y are in kpc. 



Fig. 18. Exp. F: Gray-colour maps do pattern speed Q p (r) 
(km/s/kpc) of the m=l spiral wave determined from the total 
potential and from the surface density of the gas as a function 
of the radius (in kpc) . 

The power spectrum analysis, using the amplitude and 
the phase of the Fourier components of the total potential 
and of the surface density of the gas and of the stars, is 
another method that allow us to evidence the presence of 
the spiral modes in our simulations. This analysis give us 
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the pattern speed of all present spiral perturbations and 
their relative positions. The figure [l^ displays gray-colour 
maps of the obtained angular velocity, from the poten- 
tial and from the surface density of the gas (Exp. F) as 
a function of the radius, for the m=l mode, using a such 
method. In these figure we can observe that the m=l mode 
(f2 p (r)=400km/s/kpc), identified in the distribution of the 
gas density, it is more pronounced between the centre and 
the radius 3 kpc, extending however until the edge. The 
same behaviour can be verified in this figure for the po- 
tential, indicating that although this mode be more easily 
identified in the central region of the galaxy, it produces a 
perturbation in the total potential that extends until the 
edge. 



Fig. 20. Exp. Fl: linear isodensity contours of the stars (left 
side) and gas (right side) at epochs t=20, 80, 140 and 200 (in 
10 7 yrs). X and Y are in kpc. 



Fig. 19. Exp. F: the strength of the m=l spiral perturbation 
Pi as a function of time (in 10 7 yrs). 



Figure [T^ presents the strength of the m=l (pi) mode 
as a function of time. From this, we have computed the 
growth rate of the m=l wave in the central region, assum- 
ing an exponentially growing mode. The growth time-scale 
obtained is approximately 3 id, at the radius of 3 kpc. This 
corresponds to the predictions of the SLING mechanism 
(Adams et al. 1989, Shu et al. 1990). Also, the fact that 
the m=l behaviour described here occurs only when the 
gas mass fraction is high, is compatible with this mecha- 
nism, which is essentially gaseous in its origin. 

5. Viscosity 

The dynamical role of the dissipative gas component is 
essential to account for the variety of structures observed 
in spiral galaxies. But the actual amount of dissipation 
is not well known, and worries have been expressed that 
the numerical treatment of the gas component could in- 



Fig. 21. Exp. Fl: the strength of the m=l spiral perturbation 
pi as a function of time (in 10 7 yrs). The m=l strength has 
been averaged over radii to 3 kpc. 

troduce artificial angular momentum transfer, through a 
too large viscosity. In the real interstellar medium, the 
viscosity is due to macroscopic turbulent motions, that 
can be supersonic with a lOkm/s velocity dispersion (e.g. 
Dickey et al 1990). But the corresponding viscous torques 
are quite weak at large-scale, and the time-scale for vis- 
cous transport of the galactic gaseous disk is larger than 
a Hubble time (e.g. Lin & Pringle 1987, Combes 1991). 
Gravitational torques are always dominant, and can trans- 
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Fig. 22. Exp. F and Fl : surface density (in code units) of 
stars and gas as a function of radius (kpc) at four different 
epochs t=20, 80, 140 and 200 (in 10 7 yrs). The letters at the 
bottom left corner of each frame indicate (s) stars and (g) gas. 



fer angular momentum on time-scales much shorter than 
the Hubble time. Fortunately, this makes gas-dynamical 
simulations more independent of the adopted physics of 
the interstellar medium. 

We therefore would like to check that the beam-scheme 
used here does not introduce too large viscous torques. In 
order to evaluate quantitatively the numerical viscosity of 
the method we have carried out further control experi- 
ments. The artificial viscosity must be proportional to the 
size of the grid cell where the surface density is computed. 
In the experiment Fl we have used the same initial con- 
ditions as in experiment F but the grid is now composed 
of four times more cells (512 x 512 cells) with four times 
smaller surface each. In figure |2^ we present the linear iso- 
density contours of the stars and gas for Exp. Fl at four 
different epochs. As described above for the former exper- 
iments, we calculate for experiment Fl the growth rate of 
the m=l spiral mode from the time evolution of its am- 
plitude pi plotted in figure [2l]. The corresponding growth 
time-scale is approximately 23 td, i.e. in this experiment 
the m=l spiral mode needs a much longer time-scale to 
develop, compared to Exp. F. The analysis of the surface 
density distribution of stars and g function of radius 

and at different epochs, as shown in figure ^2|, shows that 
although initially (t=20) the viscosity is larger in Exp. F, 
driving a higher gas flow towards the center, at following 
epochs (t=80) the surface density distributions of Exp. F 
and Fl evolve in a similar manner. After that, from t=100, 
their behavior diverge from one another. At t=140, the gas 
density is more concentrated in Exp. Fl, although the nu- 
merical viscosity is there lower. This must be due to the 



Fig. 23. Exp. F and Fl : radial velocity (top) and equivalent 
viscosity of the g function of time (bottom), at radii 3, 

7 and 11 kpc. 



Fig. 24. Exp. F and Fl : radial velocity (top) and equivalent 
viscosity of the gas as a function of radius (bottom), at epochs 
t=20, 80, 140 and 200 (in 10 7 yrs). 



fact that the m=l arrived later, so the gravity torques had 
more time to act than in Exp. F. 

In order to interpret this behavior, we calculated the 
equivalent "gravitational viscosity" of the gas as proposed 
by von Linden et al. (1995). From this point of view, 
the combination of gas dissipation and gravity torques 
from non-axisymmetric gravitational instabilities develop- 
ing in the disks, is considered as an equivalent "gravita- 
tional viscosity" . The gravity torques produce an angular- 
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momentum re-distribution and an evolution of the sur- 
face density of the gas. By monitoring this surface density 
one can derive the equivalent viscosity quantitatively. This 
"gravitational viscosity" does not correspond to a classi- 
cal gaseous viscosity, but it is a measure of the efficiency 
with which gravity torques transfer angular momentum 
outward and matter inward in the system. In figure ^3] 
is plotted (at top) the behavior of the radial velocity of 
the function of time at three different radii (3, 

7 and 11 kpc) for experiments F and Fl and (at bottom) 
the behavior of the corresponding "gravitational viscos- 
ity". For experiment F the viscosity has a positive value 
initially, just before the development of the m=l pertur- 
bation, and converges later on towards a zero mean value. 
In experiment Fl, the viscosity remains large for a longer 
time, which corresponds to the delay in the development 
of the m=l perturbation with respect to Exp. F. To better 
visualize the viscosity behavior in these experiments, we 
show in figure [24| the radial velocity and the gravitational 
viscosity as a function of radius at four different epochs. 
While the viscosity in Exp. F is negative at the center 
and tends to zero at larger radii after the development of 
the m=l mode, in experiment Fl the viscosity remains 
positive. 

In summary, although the numerical viscosity is 
slightly higher in Exp. F, this is largely compensated by 
the larger gravitational viscosity observed in Exp. Fl. 
When the grid cells are smaller, the gas can condense in 
cold clumps more easily (see fig |2(]), and the lumpiness 
of the gas has a significant gravitational effect. Although 
the numerical viscosity in itself does not alter significantly 
the results, the dynamic range and spatial scale of gaseous 
condensations, and therefore the physics of the ISM, plays 
a significant role. 

6. Summary and Conclusions 

We have investigated the role of gas dynamics in the 
gravitational instabilities of galaxy disks, through self- 
gravitating simulations of the coupled stellar and gaseous 
components. In the present paper, we have focussed on the 
diffuse gas component, that can be considered as a fluid, 
subject to pressure forces and shocks. 

Spiral structure is only transient in purely stellar disks. 
In the absence of a cooling mechanism the disk is heated 
by the gravitational instabilities themselves, and reaches 
a hot steady-state, that prevents the development of new 
instabilities. If a bar has developped, it can remain for a 
Hubble time. In presence of gas, even with as low a mass 
fraction as a few percents, the disk is continuously cooled 
down, and responsive to new dynamical instabilities. Due 
to its dissipative character, the gas acts as a regulator. 
If its conversion into stars and its continuous accretion 
is considered, the gas can prolong the phase of dynamical 
instability to more than a Hubble time, which can account 



for the spiral density waves observed in so many spiral 
galaxies. 

Dynamical instabilities play themselves a self- 
regulating role, since the bar gravity torques drive the 
gas towards the center, and a sufficiently massive central 
condensation begins to dissolve the bar. Also, above a cer- 
tain central gas concentration, a one armed spiral struc- 
ture can develop, which stops or slows down the radial gas 
flow and the angular momentum transfer. We have stud- 
ied the conditions required for the m=l mode appearance, 
its pattern speed and extension. The m=l perturbation is 
confined in the central part of the galaxy, and is noticeable 
only within its own outer Lindblad resonance. Its angu- 
lar speed is more than 10 times the speed of the external 
perturbations. It appears to have a purely gaseous origin, 
and corresponds to the phenomenon identified in proto- 
stellar disks by Adams et al. (1989). The mode grows in 
a few dynamical time-scales as predicted for the SLING 
mechanism of Shu et al. (1990). There is evidence from 
our simulations that stars and gas are off-centered with 
respect to the center of mass of the system, which indi- 
cates that this mechanism could be at work. 
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